{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Whence cometh the L in LU?\n",
    "\n",
    "Often, LU factorization is first presented in terms of a laborious procedure.  Getting $U$ was \"easy\", it was just Gaussian elimination.  But to get $L$, it is common to first write out the individual elimination steps as matrices, then invert them to move them to the other side, then multiply them together to get $L$.   Such a presentation makes it seem like $L$ is more trouble than it is worth, and that $A = LU$ factorization is a theoretically pretty but practically useless construction.\n",
    "\n",
    "However, nothing could be further from the truth.  It turns out that we can just \"read off\" $L$ much more simply directly from the pivot-row \"multipliers\" that we use during elimination steps, so that $L$ is an extremely useful and practical record of the elimination steps—indeed, a necessary record, in order to preserve all of the information from $A$ and to solve $Ax=b$ for new right-hand-sides $b$. \n",
    "\n",
    "To see some examples this, let's first write a Julia function to perform Gaussian elimination (without row swaps!) and print out all of the steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "print_gauss (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# perform Gaussian elimination of A without row swaps, returning U,\n",
    "# while printing a message for each elimination step.\n",
    "function print_gauss(A)\n",
    "    m = size(A,1) # number of rows\n",
    "    U = copy!(similar(A, typeof(inv(A[1,1]))), A)\n",
    "    for j = 1:m   # loop over m columns\n",
    "        for i = j+1:m   # loop over rows below the pivot row j\n",
    "            # subtract a multiple of the pivot row (j)\n",
    "            # from the current row (i) to cancel U[i,j] = Uᵢⱼ:\n",
    "            ℓᵢⱼ = U[i,j]/U[j,j]\n",
    "            println(\"subtracting $ℓᵢⱼ × (row $j) from (row $i)\")\n",
    "            U[i,:] = U[i,:] - U[j,:] * ℓᵢⱼ\n",
    "            U[i,j] = 0 # store exact zero to compensate for roundoff errors\n",
    "        end\n",
    "    end\n",
    "    return U\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's try it on a randomly chosen $5 \\times 5$ matrix $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Matrix{Int64}:\n",
       "  4  -2  -7  -4  -8\n",
       "  9  -6  -6  -1  -5\n",
       " -2  -9   3  -5   2\n",
       "  9   7  -9   5  -8\n",
       " -1   6  -3   9   6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [4  -2  -7  -4  -8\n",
    "     9  -6  -6  -1  -5\n",
    "    -2  -9   3  -5   2\n",
    "     9   7  -9   5  -8\n",
    "    -1   6  -3   9   6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "subtracting 2.25 × (row 1) from (row 2)\n",
      "subtracting -0.5 × (row 1) from (row 3)\n",
      "subtracting 2.25 × (row 1) from (row 4)\n",
      "subtracting -0.25 × (row 1) from (row 5)\n",
      "subtracting 6.666666666666667 × (row 2) from (row 3)\n",
      "subtracting -7.666666666666667 × (row 2) from (row 4)\n",
      "subtracting -3.6666666666666665 × (row 2) from (row 5)\n",
      "subtracting -1.2442748091603053 × (row 3) from (row 4)\n",
      "subtracting -0.4732824427480916 × (row 3) from (row 5)\n",
      "subtracting 33.495145631067295 × (row 4) from (row 5)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5×5 Matrix{Float64}:\n",
       " 4.0  -2.0   -7.0    -4.0        -8.0\n",
       " 0.0  -1.5    9.75    8.0        13.0\n",
       " 0.0   0.0  -65.5   -60.3333    -88.6667\n",
       " 0.0   0.0    0.0     0.262087   -0.659033\n",
       " 0.0   0.0    0.0     0.0        31.7767"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_gauss(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In comparison, here is the LU factorization of $A$ from the built-in `lu` function, with row-swaps disabled:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Matrix{Float64}:\n",
       " 4.0  -2.0   -7.0    -4.0        -8.0\n",
       " 0.0  -1.5    9.75    8.0        13.0\n",
       " 0.0   0.0  -65.5   -60.3333    -88.6667\n",
       " 0.0   0.0    0.0     0.262087   -0.659033\n",
       " 0.0   0.0    0.0     0.0        31.7767"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L, U = lu(A, NoPivot())\n",
    "U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same $U$ matrix!  Now let's look at $L$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Matrix{Float64}:\n",
       "  1.0    0.0       0.0        0.0     0.0\n",
       "  2.25   1.0       0.0        0.0     0.0\n",
       " -0.5    6.66667   1.0        0.0     0.0\n",
       "  2.25  -7.66667  -1.24427    1.0     0.0\n",
       " -0.25  -3.66667  -0.473282  33.4951  1.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the entries of $L$ below the diagonal are *exactly* the multipliers that were printed out during Gaussian elimination (the factors by which the pivot row is multiplied before it is *subtracted* from a row below it).\n",
    "\n",
    "One way to see this is to consider the matrix product $LU$, which should give $A$.  Consider, for example, the third row of this: the third row of $L$ tells us what linear combinations of the rows of $U$ gives the third row of $A$.  It says:\n",
    "\n",
    "* third row of `A` = `[-2,-9,3,-5,2]` = `-0.5` × (row 1 of U) + 6.66... × (row 2 of U) + (row 3 of U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " -2.0\n",
       " -9.0\n",
       "  3.0\n",
       " -5.0\n",
       "  2.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L[3,1] * U[1,:] + L[3,2] * U[2,:] + U[3,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this is exactly the *reverse of the elimination steps*, so of course it works.  **Putting the multipliers in $L$ is the right thing!**\n",
    "\n",
    "See also section 2.6 of the textbook for more info.\n",
    "\n",
    "Still, computing the $L$ in the LU factorization requires care to put all of the multipliers in the right place with the right sign.  It can be a pain for human beings, which is why we typically don't do it when performing Gaussian elimination by hand.  However, computers are great at this kind of tedious bookkeeping, and since keeping track of $L$ requires almost *no extra work*, computers essentially *always* figure out *both* $L$ and $U$ when doing Gaussian elimination.\n",
    "\n",
    "## Warning:\n",
    "\n",
    "It is easy to get the wrong impression here.  I will re-iterate my warning from a previous lecture: **you cannot invert an arbitrary triangular matrix just by flipping signs**. \n",
    "\n",
    "The sign-flip argument is about how elimination works, and is not about inverting an arbitrary triangular matrix.   It's not even the case that $L^{-1}$ is obtained simply by flipping signs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Matrix{Float64}:\n",
       "  1.0    0.0       0.0        0.0     0.0\n",
       "  2.25   1.0       0.0        0.0     0.0\n",
       " -0.5    6.66667   1.0        0.0     0.0\n",
       "  2.25  -7.66667  -1.24427    1.0     0.0\n",
       " -0.25  -3.66667  -0.473282  33.4951  1.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that $L^{-1}$ is *not* simply $L$ with the signs below the diagonal flipped."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another example\n",
    "\n",
    "Let's do another example, this one small enough that we can go through the calculations by hand.  We'll do Gaussian elimination on the following 3×3 invertible (non-singular) matrix:\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix} \\color{blue}{1} & 2 & 0 \\\\ 2 & 5 & 1 \\\\ -3 & 1 & -1 \\end{pmatrix} \\stackrel{r_2 - \\color{red}{2}r_1}{\\stackrel{r_3 + \\color{red}{3}r_1}{\\longrightarrow}}\n",
    "    \\begin{pmatrix} \\color{blue}{1} & 2 & 0 \\\\ 0 & \\color{blue}{1} & 1 \\\\  0 & 7 & -1 \\end{pmatrix} \\stackrel{r_3 - \\color{red}{7}r_2}{\\longrightarrow}\n",
    "    \\begin{pmatrix} \\color{blue}{1} & 2 & 0 \\\\ 0 & \\color{blue}{1} & 1 \\\\  0 & 0 & \\color{blue}{-8} \\end{pmatrix} = U\n",
    "$$\n",
    "\n",
    "(Here, \"$r_2 - 2r_1$\" denotes the elimination step \"row 2 - 2(row 1)\" etcetera.)\n",
    "\n",
    "To get $L$, we just need to write down the multipliers as we go along, with **opposite signs**, putting each multiplier in the **same column and row** as the corresponding elimination step:\n",
    "\n",
    "$$\n",
    "L = \\begin{pmatrix} 1 & & \\\\ \\color{red}{+2} & 1 & \\\\ \\color{red}{-3} & \\color{red}{+7} & 1 \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Let's check that $A = LU$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       "  1  2   0\n",
       "  2  5   1\n",
       " -3  1  -1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = [ 1  0  0 \n",
    "     +2  1  0\n",
    "     -3 +7  1 ]\n",
    "U = [ 1  2  0\n",
    "      0  1  1\n",
    "      0  0 -8 ]\n",
    "L * U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hooray, we didn't make a mistake!  We got the original $A$ back.\n",
    "\n",
    "Again, it's important to understand that what we are doing here is *not* inverting a matrix. Even for a very special matrix like $L$ (lower-triangular with 1's on the diagonal), you can't *in general* invert it simply by flipping signs below the diagonal. Let's check the inverse of $L$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       "  1.0   0.0  0.0\n",
       " -2.0   1.0  0.0\n",
       " 17.0  -7.0  1.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is *not* the same as $L$ with the signs flipped — look at the 17 in the lower-left corner.  Why is this?\n",
    "\n",
    "* Remember the formula: $A = LU$.  The matrix $L$ tells us how to get the rows of $A$ from the rows of $U$.   This uses the same multipliers as elimination, with flipped signs to reverse them, because *elimination steps always add multiples of pivot rows of **U** (not A)*!\n",
    "\n",
    "* In contrast, consider $U = L^{-1} A$.  The matrix $L^{-1}$ tells us how to get the rows of $U$ from the rows of $A$.   This is *not* the same as the elimination steps, because (except for the first row) the **pivot rows change during elimination**.   For example, in the second elimination step, we subtracted 7 times row 2 of *U* (0 1 1) from the third row of U, which is *not* the same as subtracting 7 times row 2 of *A* (2 5 1).\n",
    "\n",
    " - To get the third row of U from A above, we can *re-do* (not reverse) the elimination steps as follows: (3rd row of U) = (3rd row of A) + 3(1st row of A) - 7(2nd row of **U**) = (3rd row of A) + 3(1st row of A) - 7×((2nd row of A) - 2(1st row or A)) = (3rd row of A) + 17(1st row of A) - 7×(2nd row of A).  The coeefficients 1, 17 (= 3 + 7×2), and -7 here are exactly the third row of $L^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using LU factorizations\n",
    "\n",
    "Lots of things that you might want to do with a matrix $A$ become easier once you have the $A=LU$ factorization.  Most importantly, it becomes much easier to solve systems of equations.\n",
    "\n",
    "(Exactly *how much* easier is something we'll quantify later.  Short answer: for an $n \\times n$ matrix $A$, it takes around $n^3$ operations to perform Gaussian elimination to get $U$ and $L$, but subsequently solving for $x$ by takes only around $n^2$ operations.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving Ax=b\n",
    "\n",
    "When we do Gaussian elimination by hand, we convert $Ax = b$ to $Ux = c$ by performing the same elimination steps on $b$ to get $c$ as we performed on $A$ to get $U$.  Often, we do this by \"augmenting\" the matrix $A$ with the right-hand side $b$.  This makes it easier (*for hand calculation*) to keep track of what operations to do on $b$. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Int64}:\n",
       " -7\n",
       "  0\n",
       "  1\n",
       " -6\n",
       "  2"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [4  -2  -7  -4  -8\n",
    "     9  -6  -6  -1  -5\n",
    "    -2  -9   3  -5   2\n",
    "     9   7  -9   5  -8\n",
    "    -1   6  -3   9   6] # our random 5×5 matrix from before\n",
    "b = rand(-9:9, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       "   -7.0\n",
       "   15.75\n",
       " -107.49999999999999\n",
       "   -3.259541984732806\n",
       "  116.30097087378371"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, U_and_c = lu([A b], NoPivot()) # eliminate augmented matrix (without row swaps)\n",
    "U = UpperTriangular(U_and_c[:, 1:end-1]) # all but last column is U\n",
    "c = U_and_c[:, end] # last column is c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can solve $Ux = c$ by backsubstitution (`U \\ c`), and it should give the same answer (up to roundoff error) as `A \\ b`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×2 Matrix{Float64}:\n",
       "  2.64986    2.64986\n",
       "  1.79835    1.79835\n",
       " -0.334555  -0.334555\n",
       " -3.23373   -3.23373\n",
       "  3.65995    3.65995"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[U\\c  A\\b] # print them side by side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the computer doesn't do this: on a computer, you **almost never augment the matrix with the right-hand-side**.  Instead, you:\n",
    "\n",
    "1. Factor $A = LU$ by Gaussian elimination (not including row swaps, discussed below!), giving $Ax = b \\implies LUx = L(Ux) = b$\n",
    "2. Let $c = Ux$.  Solve $Lc = b$ for $c$ by forward-substitution.  (*Note:* this is especially easy because $L$ has only 1's on the diagonal, meaning that there are no divisions.)\n",
    "3. Solve $Ux = c$ for $x$ by backsubstitution.\n",
    "\n",
    "The key point to realize is that solving $Lc = b$ for $c$ involves *exactly the same elimination steps* as if you had augmented the matrix with $b$ during Gaussian elimination.   The bookkeeping is more tedious for a human, but computers are good at bookkeeping, and there turn out to be several practical advantages for computer software to separate solving for $LU$ and solving for $c$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       "   -7.0\n",
       "   15.75\n",
       " -107.49999999999999\n",
       "   -3.259541984732803\n",
       "  116.30097087378363"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L, U = lu(A, NoPivot()) # Gaussian elimination without row swaps\n",
    "c = L \\ b # solve Lc = b for c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same $c$ as before!\n",
    "\n",
    "Let's write a little program to write out the steps of forward-substitution so that we can see that they are indeed the elimination steps from before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c[1] = b[1] = -7.0\n",
      "c[2] = b[2]- 2.25 * c[1] = 15.75\n",
      "c[3] = b[3]- -0.5 * c[1]- 6.666666666666666 * c[2] = -107.49999999999999\n",
      "c[4] = b[4]- 2.25 * c[1]- -7.666666666666666 * c[2]- -1.2442748091603053 * c[3] = -3.259541984732806\n",
      "c[5] = b[5]- -0.25 * c[1]- -3.6666666666666665 * c[2]- -0.47328244274809156 * c[3]- 33.495145631067324 * c[4] = 116.30097087378371\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       "   -7.0\n",
       "   15.75\n",
       " -107.49999999999999\n",
       "   -3.259541984732806\n",
       "  116.30097087378371"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = similar(b, Float64)\n",
    "for i = 1:length(b)\n",
    "    print(\"c[$i] = b[$i]\")\n",
    "    c[i] = b[i]\n",
    "    for j = 1:i-1\n",
    "        print(\"- $(L[i,j]) * c[$j]\")\n",
    "        c[i] = c[i] - L[i,j] * c[j]\n",
    "    end\n",
    "    println(\" = \", c[i])\n",
    "end\n",
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Julia, `A \\ b` does this whole process for you implicitly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A smaller example\n",
    "\n",
    "Let's do another example, this one small enough to do by hand, using our 3×3 example from earlier.  Let's solve:\n",
    "\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} 1 & 2 & 0 \\\\ 2 & 5 & 1 \\\\ -3 & 1 & -1 \\end{pmatrix}}_{A = LU} \\underbrace{\\begin{pmatrix} x_1 \\\\ x_2 \\\\ x_3 \\end{pmatrix}}_x = \n",
    "\\underbrace{\\begin{pmatrix} 5 \\\\ 15 \\\\ -4 \\end{pmatrix}}_b\n",
    "$$\n",
    "\n",
    "First we solve $Lc = b$ by forward substitution:\n",
    "\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} 1 &  &  \\\\ 2 & 1 &  \\\\ -3 & 7 & 1 \\end{pmatrix}}_{L} \\underbrace{\\begin{pmatrix} c_1 \\\\ c_2 \\\\ c_3 \\end{pmatrix}}_c = \n",
    "\\underbrace{\\begin{pmatrix} 5 \\\\ 15 \\\\ -4 \\end{pmatrix}}_b\n",
    "$$\n",
    "\n",
    "This immediately gives $c_1 = 5$ from the first row, $2c_1 + c_2 = 10 + c_2 = 15 \\implies c_2 = 5$ from the second row, and $-3c_1 + 7c_2 + c_3 = -15 + 35 + c_3 = -4 \\implies c_3 = -24$  from the third row.  If you look carefully, and remember the Gaussian-elimination steps that we did on $A$, you'll see that these are in fact exactly the same elimination steps applied to $c$!\n",
    "\n",
    "Let's check this with `c = L \\ b`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       "  1  2   0\n",
       "  2  5   1\n",
       " -3  1  -1"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = [5, 15, -4]\n",
    "L = [ 1  0  0 \n",
    "     +2  1  0\n",
    "     -3 +7  1 ]\n",
    "U = [ 1  2  0\n",
    "      0  1  1\n",
    "      0  0 -8 ]\n",
    "A = L * U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       "   5.0\n",
       "   5.0\n",
       " -24.0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = L \\ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good, it matches our hand calculation!  Now we solve\n",
    "\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} 1 & 2 & 0 \\\\  & 1 & 1 \\\\  &  & -8 \\end{pmatrix}}_{U} \\underbrace{\\begin{pmatrix} x_1 \\\\ x_2 \\\\ x_3 \\end{pmatrix}}_x = \n",
    "\\underbrace{\\begin{pmatrix} 5 \\\\ 5 \\\\ -24 \\end{pmatrix}}_c\n",
    "$$\n",
    "\n",
    "by backsubstitution. The third row gives $-8 x_3 = -24 \\implies x_3 = 3$.  The second row gives $x_2 + x_3 = x_2 + 3 = 5 \\implies x_2 = 2$, and the first row gives $x_1 + 2x_2 = x_1 + 4 = 5 \\implies x_1 = 1$.  Let's check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 1.0\n",
       " 2.0\n",
       " 3.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U \\ c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, of course, this matches the solution `A \\ b` of the original system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 1.0000000000000002\n",
       " 2.0\n",
       " 2.999999999999999"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A \\ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple right-hand sides and AX = B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that we need to solve $Ax=b$ for **multiple right-hand sides** $b_1$, $b_2$, and so on.   Once we have computed $A=LU$ by Gaussian elimination, we can *re-use* $L$ and $U$ to solve each new right-hand side:\n",
    "\n",
    "1. Find $A = LU$ by Gaussian elimination\n",
    "2. Solve $Ax_1 = b_1$ by `x₁ = U \\ (L \\ b₁)`\n",
    "3. Solve $Ax_1 = b_2$ by `x₂ = U \\ (L \\ b₂)`\n",
    "4. etcetera\n",
    "\n",
    "Since solving triangular systems of equations ($L$ or $U$) is easy, this way we only do the hard/expensive part (Gaussian elimination once).\n",
    "\n",
    "Julia provides a shorthand for this process, so you don't have to worry about $L$ and $U$ and explicit forward/backsubstitution.  Instead, you compute `LU = lufact(A)`, which creates an \"LU factorization object\" `LU` that internally stores $L$ and $U$ in a compressed format (along with any permutations/row swaps as discussed below), and then you can do `LU \\ b` for each new right-hand side and it will do the (fast) triangular solves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       "  5  10   15\n",
       " 15  30   45\n",
       " -4  -8  -12"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = [b 2b 3b] # the solutions should be just x = (1,2,3), 2x, and 3x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 1.0  2.0  3.0\n",
       " 2.0  4.0  6.0\n",
       " 3.0  6.0  9.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A \\ B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 1.0  2.0  3.0\n",
       " 2.0  4.0  6.0\n",
       " 3.0  6.0  9.0"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A^-1 * B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearAlgebra.LU{Float64, Matrix{Float64}}\n",
       "L factor:\n",
       "3×3 Matrix{Float64}:\n",
       "  1.0       0.0       0.0\n",
       " -0.666667  1.0       0.0\n",
       " -0.333333  0.411765  1.0\n",
       "U factor:\n",
       "3×3 Matrix{Float64}:\n",
       " -3.0  1.0      -1.0\n",
       "  0.0  5.66667   0.333333\n",
       "  0.0  0.0      -0.470588"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LU = lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Matrix{Float64}:\n",
       " 1.0  1.0\n",
       " 2.0  2.0\n",
       " 3.0  3.0"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[LU\\b  A\\b] # print them side by side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equivalently, if we let $B = (b_1 \\; b_2 \\; \\cdots)$ be the matrix whose columns are the right-hand sides, and $X = (x_1 \\; x_2 \\; \\cdots)$ be the matrix whose columns are the solutions, then solving $Ax_1 = b_1$, $Ax_2 = b_2$, … is equivalent to solving $AX = B$, because $AX = (Ax_1 \\; Ax_2 \\; \\cdots)$ in the \"matrix × columns\" picture of matrix multiplication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 2.0000000000000004\n",
       " 4.0\n",
       " 5.999999999999998"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LU \\ 2b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 1.0  2.0  3.0\n",
       " 2.0  4.0  6.0\n",
       " 3.0  6.0  9.0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LU \\ B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It gives the same answer!  On a computer, solving for a bunch of right-hand sides at once by `A \\ B` is often **more efficient** than solving them one by one (for technical reasons involving the speed of memory access).   Conceptually, it is often convenient to think of many right-hand sides and solutions together, in a matrix, rather than separately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Row swaps vs. LU\n",
    "\n",
    "Occasionally, we may encounter a zero in the pivot position.  Sometimes this means that the equations are **singular** (may have no solutions) — we will talk more about this later.  However, as long as there is a nonzero value *below* the pivot, we can fix the problem by **swapping rows** (which just corresponds to re-ordering the equations).\n",
    "\n",
    "For example:\n",
    "\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{rrr|r}\n",
    "\\boxed{1} & 3 & 1 & 9 \\\\\n",
    "1 & 3 & -1 & 1 \\\\\n",
    "3 & 11 & 6 & 35\n",
    "\\end{array}\\right]\n",
    "\\stackrel{r_2 - \\color{red}{1}r_1}{\\stackrel{r_3 + \\color{red}{3}r_1}{\\longrightarrow}}\n",
    "\\left[\\begin{array}{rrr|r}\n",
    "\\boxed{1} & 3 & 1 & 9 \\\\\n",
    "0 & 0 & -2 & -8 \\\\\n",
    "0 & 2 & 3 & 8\n",
    "\\end{array}\\right]\n",
    "\\stackrel{\\mbox{swap } r_3 \\leftrightarrow r_2}{\\longrightarrow}\n",
    "\\left[\\begin{array}{rrr|r}\n",
    "\\boxed{1} & 3 & 1 & 9 \\\\\n",
    "0 & \\boxed{2} & 3 & 8 \\\\\n",
    "0 & 0 & \\boxed{-2} & -8\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "where in the second step we swapped the second and third rows to get a nonzero pivot in the second row.\n",
    "\n",
    "At this point we can again solve bottom-up by backsubstitution:\n",
    "\n",
    "$$\n",
    "-2x_3 = 8 \\implies x_3 = 4 \\\\\n",
    "2x_2 + 3x_3 = 8 = 2x_2 + 12 \\implies x_2 = -2 \\\\\n",
    "x_1 + 3x_2 + x_3 = 9 = x_1 -6 + 4 \\implies x_3 = 11\n",
    "$$\n",
    "\n",
    "Of course, the computer can get the answer much more quickly and easily:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 1   3   1\n",
       " 1   3  -1\n",
       " 3  11   6"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [ 1 3  1\n",
    "      1 3 -1\n",
    "      3 11 6 ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Int64}:\n",
       "  9\n",
       "  1\n",
       " 35"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = [9, 1, 35]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 11.000000000000005\n",
       " -2.0000000000000013\n",
       "  4.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A \\ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But what does this mean for LU factorization?  If we write out the elimination steps (laboriously) as matrix multiplications, then the second step is expressed by a [permutation matrix](https://en.wikipedia.org/wiki/Permutation_matrix) $P_2$ that swaps the 2nd and third rows:\n",
    "\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} 1 & 3 & 1 \\\\ & 2 & 3 \\\\ & & -2 \\end{pmatrix}}_U =\n",
    "\\underbrace{\\begin{pmatrix} 1 &  &  \\\\ &  & 1 \\\\ & 1 &  \\end{pmatrix}}_{P_2}\n",
    "\\underbrace{\\begin{pmatrix} 1 &  &  \\\\ -1 & 1 &  \\\\ -3 &  & 1  \\end{pmatrix}}_{E_1}\n",
    "\\underbrace{\\begin{pmatrix} 1 & 3 & 1 \\\\ 1 & 3 & -1 \\\\ 3 & 11 & 6 \\end{pmatrix}}_A\n",
    "$$\n",
    "\n",
    "But if we multiply $P_2 E_1$, the result is no longer lower triangular:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       "  1  0  0\n",
       " -3  0  1\n",
       " -1  1  0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P2 = [1 0 0\n",
    "      0 0 1\n",
    "      0 1 0]\n",
    "E1 = [1 0 0\n",
    "     -1 1 0\n",
    "     -3 0 1]\n",
    "P2 * E1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason for this is straightforward: the lower-triangular structure of the elimination matrix $E$ (and hence $L = E^{-1}$) came from the fact that we **always subtracted earlier rows from later rows** during elimination, and this **is no longer true** if we did row swaps.\n",
    "\n",
    "Does this mean that the whole concept of LU factorization goes out the window?  Is $A$ no longer a product of \"nice\" triangular matrices?  **No.**  There are two tricks.\n",
    "\n",
    "**Trick #1**: if we had re-ordered the rows of $A$ *in the beginning*, then we wouldn't have needed row swaps:\n",
    "\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} 1 &  &  \\\\ &  & 1 \\\\ & 1 &  \\end{pmatrix}}_{P}\n",
    "\\underbrace{\\begin{pmatrix} 1 & 3 & 1 \\\\ 1 & 3 & -1 \\\\ 3 & 11 & 6 \\end{pmatrix}}_A\n",
    "=\n",
    "\\underbrace{\\begin{pmatrix} 1 & 3 & 1 \\\\ 3 & 11 & 6 \\\\ 1 & 3 & -1  \\end{pmatrix}}_{PA}\n",
    "=\n",
    "\\underbrace{\\begin{pmatrix} 1 &  &  \\\\ 3 & 1 &  \\\\ 1 & & 1 \\end{pmatrix}}_L\n",
    "\\underbrace{\\begin{pmatrix} 1 & 3 & 1 \\\\ & 2 & 3 \\\\ & & -2 \\end{pmatrix}}_U \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 1   3   1\n",
       " 3  11   6\n",
       " 1   3  -1"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = [1 0 0\n",
    "     3 1 0\n",
    "     1 0 1]\n",
    "U = [1 3 1\n",
    "     0 2 3\n",
    "     0 0 -2]\n",
    "L*U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is called the \"PA=LU\" factorization, and is what *real* numerical linear-algebra software actually computes, as we discuss more below.\n",
    "\n",
    "**Trick #2:** How do you know what re-ordering of A to use *in advance?*  It turns out that you don't need to do it in advance — there is a tricky book-keeping technique to \"work backwards\" to figure out what P and L as you go along.  Whenever you do a row swap, there is a way to \"commute\" it through the elimination steps to figure out what L *would* have been if you had *started* with that permutation!\n",
    "\n",
    "I *won't* go over this book-keeping trick in 18.06!   It's more a technique for computer implementations, and is described in many textbooks on numerical linear algebra.   I will **never** ask you to compute the \"PA=LU\" factorization by hand.   But it is important to know that:\n",
    "\n",
    "* PA=LU is the factorization that you get in *practice* from *real* Gaussian elimination steps.\n",
    "* U is the output of elimination, as usual, and L and P are just a record of the elimination and row-swap steps (just book-keeping, no arithmetic required).\n",
    "* If you are *given* PA=LU, we need to know how to *use* it, e.g. to solve Ax=b."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Row swaps and PA = LU in practice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up to now, we have mostly ignored the possibility of row swaps.  Row swaps may be *required* if you encounter a zero pivot (assuming there is a nonzero value below it in the same column), but this is often unlikely to occur in practice (especially for random matrices!).\n",
    "\n",
    "However, even as in our examples where no row swaps were *required*, a computer will often do them *anyway*, in order to minimize roundoff errors.  It turns out that roundoff errors (the computer only keeps about 15–16 significant digits) can be disastrous if the pivot is merely very small.  So, the computer swaps rows to **make the pivot as big as possible**, as strategy called *partial pivoting*.  As a result, the `lu` function in Julia returns *three* things: $L$, $U$, and the permutation $p$ giving the **re-ordering of the rows** of $A$ that is needed.  For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       "  1  2   0\n",
       "  2  5   1\n",
       " -3  1  -1"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# matrix from earlier example, does not *require* row swaps\n",
    "\n",
    "A = [ 1  2   0\n",
    "      2  5   1\n",
    "     -3  1  -1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       "  1.0       0.0       0.0\n",
       " -0.666667  1.0       0.0\n",
       " -0.333333  0.411765  1.0"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L, U, p = lu(A)\n",
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " -3.0  1.0      -1.0\n",
       "  0.0  5.66667   0.333333\n",
       "  0.0  0.0      -0.470588"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that this is **not the same as what we got by hand!**  Julia (really, the underlying [LAPACK linear-algebra library](https://en.wikipedia.org/wiki/LAPACK)) re-orders the rows *even though it doesn't have to*.  (If you want to look it up, the algorithm is called \"**partial pivoting**\".)\n",
    "\n",
    "To know what re-ordering Julia's `lu` function did, we can look at the `p` vector returned by `lu(A)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Int64}:\n",
       " 3\n",
       " 2\n",
       " 1"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`p` says tells you in what order we should put the rows of $A$ to match the product $LU$: we should re-order `A` to put row 3 first, then row 2, then row 1.  We can do this in Julia easily by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " -3  1  -1\n",
       "  2  5   1\n",
       "  1  2   0"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[p,:] # A with the rows in order p = PA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should match $LU$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " -3.0  1.0  -1.0\n",
       "  2.0  5.0   1.0\n",
       "  1.0  2.0   0.0"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L*U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The computer only stores a list of numbers for `p` because that is the most efficient way to store and work with the permutation.  However, for algebraic manipulations it is often convenient to think of this as a **permutation matrix** $P$ multiplying $A$.  Since $P$ re-orders the *rows* of $A$, it must multiply $A$ on the *left*.  Constructing $P$ is easy: we just **apply the same row permutation to the identity matrix I**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "permutation_matrix (generic function with 1 method)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# construct a permutation matrix P from the permutation vector p\n",
    "permutation_matrix(p) = I(length(p))[p,:] # just reorder the rows of I (returned by I(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 SparseArrays.SparseMatrixCSC{Bool, Int64} with 3 stored entries:\n",
       " ⋅  ⋅  1\n",
       " ⋅  1  ⋅\n",
       " 1  ⋅  ⋅"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = permutation_matrix(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " -3  1  -1\n",
       "  2  5   1\n",
       "  1  2   0"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, LU factorization with row swaps corresponds to a factorization\n",
    "\n",
    "$$\n",
    "PA = LU \\Longleftrightarrow A = P^{-1} LU \\Longleftrightarrow A^{-1} = U^{-1} L^{-1} P\n",
    "$$\n",
    "\n",
    "Now, to solve $Ax = b$, a more complete process is:\n",
    "\n",
    "1. Factor $PA = LU$\n",
    "2. Compute $x =  U^{-1} L^{-1} P b$, multiplying left to right:\n",
    "   1. First, multiply $Pb$: permute $b$ according to $P$, i.e. compute `b[p]` in Julia\n",
    "   2. Second, multiply $c = L^{-1} (P b)$, i.e. solve $Lc = Pb$ for $c$ by forward-substitution\n",
    "   3. Third, multiply $x = U^{-1} c$, i.e. solve $Ux = c$ for $x$ by backsubstitution.\n",
    "\n",
    "Of course, Julia does all of this for you automatically with `A \\ b` or `lufact(A) \\ b`, but we can do it manually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       "  -2.25\n",
       "   5.625\n",
       " -22.625"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A \\ b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Int64}:\n",
       " 35\n",
       "  1\n",
       "  9"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b[p] # permute b (equivalent to Pb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       "  -2.25\n",
       "   5.625\n",
       " -22.625"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = L \\ b[p] # solve Lc = Pb = b[p]\n",
    "x = U \\ c # solve Ux = c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearAlgebra.LU{Float64, Matrix{Float64}}\n",
       "L factor:\n",
       "3×3 Matrix{Float64}:\n",
       "  1.0       0.0       0.0\n",
       " -0.666667  1.0       0.0\n",
       " -0.333333  0.411765  1.0\n",
       "U factor:\n",
       "3×3 Matrix{Float64}:\n",
       " -3.0  1.0      -1.0\n",
       "  0.0  5.66667   0.333333\n",
       "  0.0  0.0      -0.470588"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LU = lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       "  -2.25\n",
       "   5.625\n",
       " -22.625"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LU \\ b # applies permutation, then forward-sub for L, then back-sub for U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hooray, this is the same answer as above!\n",
    "\n",
    "One final point often confuses people here, if you think carefully about the above process.   By writing $PA = LU$, it seems like you *first* decide on the row-reordering of $A$, and *then* compute the LU factorization of $PA$.  But how do you know the proper row-reordering *before* you do elimination?  In fact, this is an illusion: the computer figures out the row-reordering as it goes along (partial pivoting as described above), but it then cleverly works backwards to figure out what reordering it *should* have done in the beginning!"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.7.1",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
